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Abstract. Power series representations for special functions are computationally 
satisfactory only in the vicinity of the expansion point. Thus, it is an obvious idea to 
use instead Pade approximants or other rational functions constructed from sequence 
transformations. However, neither Pade approximants nor sequence transformation 
utilize the information which is available in the case of a special function - all 
power series coefficients as well as the truncation errors are explicitly known - in 
an optimal way. Thus, alternative rational approximants, which can profit from 
additional information of that kind, would be desirable. It is shown that in this 
way a rational approximant for the digamma function can be constructed which 
possesses a transformation error given by an explicitly known series expansion. 
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1. Introduction 

Power series are among the most important tools of calculus. For ex- 
ample, they are extremely useful for the construction of solutions to 
differential equations. Accordingly, many special functions are defined 
and computed via power series. 

However, from a purely numerical point of view, a power series 
representation is a mixed blessing. Power series converge well only in 
the vicinity of the expansion point. Further away, they converge slowly 
or even diverge. Consequently, the defining power series alone normally 
does not suffice for an efficient and reliable computation of a special 
function. 

In applied mathematics and in theoretical physics, Pade approxi- 
mants have become the standard tool to overcome convergence prob- 
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lems with slowly convergent or divergent power series [2J. Therefore, it 
looks like an obvious idea to use them for the computation of special 
functions. 

Pade approximants are defined as solutions of a system of nonlinear 
equations [2], although they are in practice more often computed by 
recursive algorithms, for example by Wynn's epsilon algorithm All 
these algorithms only need the input of the numerical values of the lead- 
ing series coefficients. No further information about the function, which 
is to be approximated, is needed. This is a very advantageous feature, 
in particular if apart from a finite number of series coefficients very 
little else is known, and it has undoubtedly contributed significantly to 
the popularity of Pade approximants and their practical usefulness. 

If, however, we want to compute a special function, we are in a much 
better situation. Not only do we know explicitly all coefficients of the 
power series, but we are also able to write down at least formally explicit 
expressions for the truncation errors. An approximation scheme for a 
special function should be able to benefit from additional information 
of that kind, but Pade approximants - due to their very nature - 
cannot. Thus, valuable information is wasted, and Pade approximants 
are in the case of special functions normally less effective than other 
sequence transformations which can utilize information of that kind. 
For example, it was shown in ^1 El El EE] that Levin's sequence 
transformation jS] and some generalizations |121 Sections 7-9] can be 
much more effective than Pade approximants, in particular if factorially 
divergent asymptotic series for special functions have to be summed. 

The power of Levin-type transformations, which were recently re- 
viewed by Homeier (I] , is due to the fact that they use as input data not 
only the elements of the sequence to be transformed but also explicit 
truncation error estimates. In the majority of all applications, only 
some very simple truncation error estimates introduced by Levin jS] 
and Smith and Ford are used. In the case of special functions, 
however, it may well be possible to derive more sophisticated trun- 
cation error estimates which should ultimately lead to more effective 
approximations. Further research into this direction would be highly 
desirable. 

In the case of special functions, it is also possible to pursue a more 
direct approach. As discussed for example in a sequence trans- 

formation is a map, which transforms a slowly convergent or divergent 
sequence {s n }^ =0 , whose elements may for instance be the partial sums 
of an infinite series, into another sequence {s^I^Lq with hopefully bet- 
ter numerical properties. Concerning the input sequence it is assumed 
that its elements can for all n € No be partitioned into a (generalized) 
limit s and a remainder r n according to s n = s + r n . A sequence 
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transformation tries to determine and eliminate the remainders r n from 
the sequence elements s n . Unfortunately, a complete elimination of the 
remainders can normally be accomplished only in the case of more or 
less artificial model problems. Thus, the elements of the transformed 
sequence can also be partitioned according to s' n = s + r' n into the same 
(generalized) limit s and a transformed remainder r' n which is normally 
nonzero for all finite values of n. The transformation process was suc- 
cessful if the transformed remainders {r^}^ have better numerical 
properties than the original remainders {r n }^L - 

Normally, only relatively little is known about the remainders r n . 
In the case of special functions, however, the situation is much bet- 
ter: All coefficients of the power series are explicitly known, and the 
truncation errors of the partial sums of the power series are at least in 
principle also explicitly known. Thus, it should be possible to optimize 
the determination and elimination of the remainders - or equivalently 
the transformation process - by utilizing the available information as 
effectively as possible. 

It is the intention of this article to show that these goals can be 
accomplished in the case of the psi or digamma function ^ Eq. (6.3.1)] 

which is a meromorphic function with poles at z = 0, —1, —2, .... Our 
starting point is the power series representation ^ Eq. (6.3.14)] 

${l + z) = - 7 + zZ(z), (1.2a) 

oo 

Z(z) = £ C(* + 2) (-*)" , (L2b) 

which converges for \z\ < 1. Here, 7 is Euler's constant ^ Eq. (6.1.3)] 
and Q{y + 2) is a Riemann zeta function ^ Eq. (23.2.1)]. 



2. The transformation of the power series 

In this article, an explicit rational approximant for Z{z) will be con- 
structed. We only have to consider < z < 1. For z < 0, we can use the 
reflection formula — z) = ip(z) + it coth(nz) ^ Eq. (6.3.5)], and for 
z > 1, we can use the recurrence formula ip(z + 1) = ip(z) + l/z \Q Eq. 
(6.3.5)]. If the argument z is very large, the digamma function should 
of course be computed via its asymptotic expansion ^ Eq. (6.3.18)]. 
For our purposes, it is convenient to rewrite (|1.2b|) as follows: 

Z(z) = Z n {z) + K n {z) , (2.1a) 
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z n {z) = e c(^+2)(-^r, (2.1b) 

oo 

K n {z) = {-z) n+l + , + (-*)" • (2.1c) 

!/=0 

As discussed in the previous section, a rational approximant to Z{z) 
can only improve convergence if the truncation errors lZ n {z) are trans- 
formed into other truncation errors with better numerical properties. 
Thus, we first have to rewrite 7Z n (z) in such a way that we better under- 
stand its nature. This can be achieved by replacing the zeta functions 
C(n + v + 3) in (|2~Tc)) by their Dirichlet series [T[ Eq. (23.2.1)] and by 
interchanging the order of summations. The resulting inner series is a 
geometric series and can be expressed in closed form. Thus, we obtain 

ZJz) = Z(z) - (-l) n+1 T [z/(rn + l)] n+1 
U{ 1 {) { ' if (m + l)(m + z + l) 1 ' 

This relationship, which holds for all z ^ —1, —2, —3, . . ., shows that the 
partial sums Z n {z) are a special case of the following class of sequences 
with qj = z/j and cj = + z)\: 

oo 

s n = s + (-1)" +1 £ . n e N o • (2.3) 

Concerning the gj's, we assume that they all have the same sign and 
are ordered in magnitude according to 

1 > M > |tts| > • • • > |9(| > |«+i| > • • • > , (2.4) 

whereas the Cj's are unspecified coefficients. 

Wynn |18j showed that the convergence of a sequence of the type of 
(|2.3[) can be accelerated with the help of his epsilon algorithm 17 : 

e { -l = 0, = s n , nGN , (2.5a) 

4i = 4-i 1} + Vt4 n+1) - 4 n) ] , MeNo. (2.5b) 

The epsilon algorithm requires as input data only the numerical values 
of the elements of the sequence (|2.3|) . but not the values of the q/s. 
Wynn also derived asymptotic estimates for the transformation errors 
s — e& |18l Theorems 16 and 17], which were later extended by Sidi 

Although the epsilon algorithm is a very powerful accelerator for 
sequences of type of (|2.3[) - numerical studies showed that the epsilon 
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algorithm accelerates the convergence of the power series in 
much more effectively than for example Levin's transformation [3] or 
some generalizations |121 Sections 7 - 9] - it nevertheless cannot profit 
from the fact that in the case of (|2.2j) the qj are explicitly known. 
Moreover, no explicit expressions for the rational approximants or the 
transformation errors are known. Thus, we use instead the sequence 
transformation 

4 n) = T^\s ni ...,s n+k ) = fl (2.6) 

K=l L ^1k 

as the starting point for the construction of an explicit rational ap- 
proximant to Z{z). Here, E is the shift operator defined by Ef(n) = 
/(n+1). 

The sequence transformation can also be computed recursively: 



T (n) = s n , n E N , (2.7a) 

(n) = h + Qj±lh k n e Nq (2 7b) 

+ 1 + Qk+l 

Essentially identical sequence transformations were discussed by Matos 

mini 



An explicit expression for can be derived with the help of the 



elementary symmetric polynomials in n variables x\, . . ., x n , which 
are defined by the generating function 6, p. 13] 

n n 

na+^t) = E4 n) ^- (2.8) 

The substitution s = 1/t yields the equivalent generating function 





Comparison with ()2.6j) shows that possesses an explicit expression 

(k) 

involving the elementary symmetric polynomials e K in the k variables 
q K with 1 < k < k: 

k k 

4-«(9l> •••?*) E " s « H e k-is(?l> ■■■Ik) Sn+K 
k k k 

K=0 K=0 

(2.10) 
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In most practical applications, the sequence transformation is 
not particularly useful since the values of the qj's have to be explicitly 
known. If, however, this is the case and if the input data are the ele- 
ments of the sequence (|2.3j) , then it can be shown by complete induction 
in k that 

oo k 

if) = a + (-ir 1 £ c 3 n ^-f ^r +1 ■ (2.H) 

j=k+l K=l Qk + 1 

Thus, the first k exponential terms Cj(qj) n+1 in (|2.3j) are eliminated. 
Since the qj's are by assumption ordered in magnitude according to 
(|2.4() . this leads to an acceleration of convergence. Moreover, the appli- 
cation of to the elements of the sequence ()2,3j) leads for sufficiently 
large values of k to a convergent sequence if the original sequence 
diverges because the leading qj's in (|2.3j) satisfy \qj\ > 1. 

Combination of (|2.2j) and (|2.11j) yields the following explicit rational 
approximant to Z(z): 

T^(Z n (z),...,Z n+k (z)) = n l±MA Zn{z) = Z {z) 

K=l 1 +{ Z / K ) 

_(_!)«+! i^L y ^±}h ( 212) 

(z + l)k (k + m+ l) n + k + 2 (k + m + z + lf ' ' 

Here, (z + 1)^ and (m + l)fc are Pochhammer symbols. It is a remark- 
able feature of this rational approximant that its transformation error 
possesses an explicit series expansion. This is quite uncommon in the 
theory of rational approximants. Moreover, the first k poles of Z{z) 

are reproduced by , whereas the remaining 
poles at z = —k — l,—k — 2,... are reproduced by the infinite series for 
the transformation error. 

The prefactor of the infinite series in (|2.12|) can be expressed as a 
beta function B(x,y) = T (x)T (y) /T [x + y) Eq. (6.2.2)] according to 



n+fc+l z n+k+2 



{z + l) k kl 



B(z,k + 1). (2.13) 



It is possible to derive an alternative expression for the infinite series 
in (|2.12|) which more closely resembles the infinite series in (|2.2|) . from 
which it was derived. For that purpose, we write the Pochhammer 
symbol in the infinite series in ()2.12j) as a product according to (m + 
l)fc = l\ k K =i{[k + m + l] + [K-k-l]). Comparison with (|2~U|) shows that 
this is the generating function of the elementary symmetric polynomials 
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in the k variables x R = k — k — 1 with 1 < k < k. Thus, we obtain 

k k 

s(fc) 



(m + l) fe = J2 €- K ( k + m + 1 T = Y,^ ] {k + m + l) k - K . (2.14) 
Inserting this into ()2.12j) yields: 

T r [ k Hz n (z),...,Z n+k (z)) = Z{z)- (-ir +1 



(2 + 1)* 
k oo -. 

^ (fc + m + l)™+ K + 2 (A; + m + z + l) ' 1 ' 

If we now do a Taylor expansion of l/(& + m + z + 1) and introduce 
the generalized (Hurwitz) zeta function ((z, a) = J2^=o( n + o)~ z with 
-1, -2, . . . p. 22], we obtain 



E 



m=0 

oo 



{k + m + l) n+ti + 2 (k + m + z + 1) 

oo 

^ CO + + k + 3, k + 1) (-z) m (2.16) 



m=0 



and 



r( fc )(^ n (z),...,z n+fc (z)) = z[z)- (-i) n+1 



oo A; 

x E (" z ) m E g l fe) C(n + m + k + 3, jb + 1) . (2.17) 

m=0 k=0 

Further modifications of (|2.17|) are possible. For example, ordinary 
zeta functions can be introduced instead of the generalized (Hurwitz) 
zeta functions according to 

k 

C(n+m+n+3,k+l) = C(n+m+K+3) - ^(i/+l)~ n ~ m ~ K ~ 3 . (2.18) 

For larger values of k, the inner sum in Q2.17J) is likely to become nu- 
merically unstable since the k variables x K = k — k — 1 of the elementary 
symmetric polynomials el are all negative. This implies that the 
alternate in sign with increasing k. Thus, sums of the type Y^,k=o e« /« 
seem to have similar properties as sums of the type J2k=o(~ 1) k (k)/k 
which are known to be numerical unstable for larger values of k if all 
f K have the same sign. 
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3. Numerical examples 

We now want to show that the new explicit rational approximant to 
Z(z) is indeed a numerically useful tool. Thus, we apply both 
defined by (|2,6|) as well as Wynn's epsilon algorithm ()2.5j) to the partial 
sums Z n (z) defined by (EHbl) . 

(k) 

The transforms T„ were computed with the help of the recurrence 
formula (|2.7|) which should be the most effective approach. For the re- 
cursive calculation, two one-dimensional arrays t and q suffice (compare 
[T3 Sections 4.3 and 7.5]): 

t[0] «- s , (3.1a) 

t[m] <- s m , q[m] <- q m , m>l, (3.1b) 

t\m — j + 1] + q\j]t\m — j] 
t[m-j] <- i 1 -, l<j<m. 3.1c 

For each m > 0, t[0] = is used to approximate the limit of the 
input sequence. 

The argument z = 1 considered in Table El lies on the boundary 
of the circle of convergence of the power series (|1.2b)l for Z(z). The 
approximants in the last column of Table 01 were chosen according to 
[fa Eq. 4.3-6]. 

All calculations in Table 01 were done in MapleV Release 5.1 with 
an accuracy of 32 decimal digits. When the accuracy was reduced to 
16 digits, at most the last digit printed differed. Thus, the computa- 
tion of the rational approximants in Table 01 is apparently numerically 
remarkably stable. 

The results in Table 01 show that the new rational approximant and 
Wynn's epsilon algorithm produce results of virtually identical quality. 
This is also observed in the case of complex arguments. In Table 01 we 
consider z = [1 + y/3 i] /2 which again lies on the boundary of the circle 
of convergence of the power series (jl.2bj) for Z(z). 

In the case of the epsilon algorithm, we obtain in the case of z = 
[l + V3i]/2: 

- 7 + zef} = 0.285 073 441 270 305 + 0.691 215 820 928 757^.2) 
- 7 + ze W = 0.285 073 441 270 304 + 0.691 215 820 928 755(i3.3) 

All calculations in Table 01 were again done with an accuracy of 32 
decimal digits. When we reduced the accuracy to 16 digits, we observed 
as in Table 01 that at most the last digit printed differed. 

Wynn's epsilon algorithm is - as already remarked - a very powerful 
accelerator for sequences of the type of (|2.3j) . Thus, the numerical 
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Table I. Convergence of the new rational approximant for ip(l + z) with z = 1. 



n 


-7-f 


- zZ n {z) 






-7 + zJ 


(0) 
n 






-J + zel 


n-2[[n/21) 
>4n/2] 





1 


.067 


718 


1 


.067 


718 


401 


946 


694 


1 


.067 


718 


401 


946 


694 


1 


— 


.134 


339 





.466 


689 


950 


366 


896 


— 0. 


.134 


338 


501 


212 


901 


2 





.947 


985 





.426 


778 


727 


217 


411 


0. 


.435 


187 


600 


653 


266 


3 


— 


.088 


943 





.423 


160 


888 


178 


296 


0. 


.418 


415 


084 


082 


869 


4 





.928 


400 





.422 


818 


740 


326 


191 


0. 


.422 


960 


666 


980 


241 


5 


-0 


.079 


949 





.422 


787 


312 


867 


790 


0. 


.422 


747 


356 


295 


448 


6 





.924 


128 





.422 


784 


577 


276 


038 


0. 


.422 


786 


030 


269 


854 


7 


-0 


.077 


880 





.422 


784 


353 


568 


296 


0. 


.422 


784 


084 


294 


859 


8 





.923 


114 





.422 


784 


336 


420 


153 


0. 


.422 


784 


346 


626 


811 


9 


-0 


.077 


380 





.422 


784 


335 


187 


365 


0. 


.422 


784 


333 


783 


337 


10 





.922 


866 





.422 


784 


335 


104 


100 


0. 


.422 


784 


335 


156 


547 


11 


-0 


.077 


256 





.422 


784 


335 


098 


804 


0. 


.422 


784 


335 


093 


078 


12 





.922 


805 





.422 


784 


335 


098 


486 


0. 


.422 


784 


335 


098 


692 


13 


-0 


.077 


226 





.422 


784 


335 


098 


468 


0. 


.422 


784 


335 


098 


450 


14 





.922 


789 





.422 


784 


335 


098 


467 


0. 


.422 


784 


335 


098 


468 



ip(l + z) 0.422 784 335 098 467 0.422 784 335 098 467 



results presented in Tables 01 and |3 indicate that the new rational 
approximant to Z(z) is indeed numerically useful for the com- 
putation of the digamma function. Nevertheless, further improvements 
should be possible. For example, we so far completely ignored that we 
have an explicit series expansion for the transformation error accord- 
ing to ()2.12j) . The convergence of this series can also be accelerated, 
for instance by Levin's u transformation jS], but unfortunately, its 
convergence cannot be accelerated as effectively as the convergence of 
the series (jl.2b|l for Z(z). So, most desirable would be an asymptotic 
expansion for the transformation error in (|2. 12|) as k becomes large. 
However, this is beyond the scope of the present article. 
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Table II. Convergence of the new rational approximant for + z) with 

z = [l + V3 i] /2. 



n -7 + z Zn (z) -7 + zTA 









245 


+ 


1 


425 i 





245 


251 


368 


522 


580 


+ 


1 


424 


554 


689 


141 


014 


1 





846 


+ 





384 i 





245 


251 


368 


522 


580 


+ 





730 


546 


812 


820 


574 


2 


-0 


236 


+ 





384 i 





279 


460 


988 


364 


996 


+ 





691 


044 


946 


370 


787 


3 


() 


282 




1 


282 i 





284 


708 


842 


795 


361 







690 


769 


505 


446 


420 


4 





791 


+ 





401 i 





285 


084 


829 


446 


025 


+ 





691 


160 


242 


235 


571 


5 


-0 


217 


+ 





401 i 





285 


078 


145 


123 


076 


+ 





691 


213 


499 


135 


601 


6 





285 


+ 


1 


270 i 





285 


073 


844 


960 


382 


+ 





691 


216 


025 


583 


709 


7 





786 


+ 





402 i 





285 


073 


447 


077 


753 


+ 





691 


215 


856 


873 


046 


8 


-0 


215 


+ 
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285 


073 


439 


323 


115 


+ 
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215 


822 


849 


613 
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285 


+ 


1 


269 i 





285 


073 


441 


081 


160 


+ 





691 


215 


820 


893 


795 


10 





785 


+ 





403 i 





285 
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441 


265 


135 


+ 





691 


215 


820 


917 


156 


11 


-0 


215 


+ 
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285 
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441 


270 


720 
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691 


215 


820 


928 


085 


12 





285 


+ 


1 


269 i 





285 
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441 


270 


350 


+ 





691 


215 


820 


928 


754 


13 





785 


+ 





403 i 





285 


073 


441 


270 


305 


+ 





691 


215 


820 


928 


757 


14 


-0 


215 


+ 





403 i 





285 


073 


441 


270 


303 


+ 





691 


215 


820 


928 


756 


15 





285 


+ 


1 


269 i 





285 


073 


441 


270 


304 


+ 





691 


215 


820 


928 


755 



tp(l + z) 0.285 073 441 270 304 + 0.691 215 820 928 756 
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